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Abstract 

Rapid trajectory generation is crucial to prompt global warfare. To meet the 
USAF’s objective of Persistent and Responsive Precision Engagement, a rapid mission 
planning tool is required. This research creates the framework for the mission planning 
tool and provides a sample optimal trajectory which is solved using the GPOPS software 
package. GPOPS employs a Gaussian pseudospectral method to solve the non-linear 
equations of motion with both end conditions and path constraints. By simultaneously 
solving the entire trajectory based on an initial guess and small number of nodes, this 
method is ideal for generating rapid solutions. The sample case is a multi-phase 
minimum time, optimal control problem which is used to validate the planning tool. The 
developed framework includes different atmospheric models, gravity models, inclusion 
of no-flyzones and waypoints, and the ability to create a library of sample cases. This 
versatile tool can be used for either trajectory generation or mission analysis. 

The results of this research show the complexities in solving an optimal control 
problem with states that change from one phase of the problem to another. At the 
conclusions of this research multiple phases were successfully connected and solved as a 
single optimal control problem. However, the entire trajectory solution from launch to 
impact solved simultaneously, is still an objective yet to be demonstrated. The results 
found should be a solid foundation for a future mission planning tool. 
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SIMULATION AND APPLICATION OF GPOPS FOR TRAJECTORY 
OPTIMIZATION AND MISSION PLANNING TOOL 


I. Introduction 

Hollywood has romanticized the concept of the strategic super computer and war. 
Movies imply that the all knowing machines of artificial intelligence have their own 
sense of awareness and eventually can take over the world. In fact, computers are more 
capable than humans of sorting through vast amount of information in performing 
repetitive calculations. With the current advances in computing power, the ability to 
rapidly conduct complex strategic and tactical war plans, like in Hollywood storylines, is 
only a question of when. 

1.1 Motivation 

A mission planning tool is very applicable to today’s Air Force mission. A tool, 
which can rapidly produce optimal trajectories has both operational and research 
potential. Operationally, this tool can be put to use in support the new Global Strike 
Command mission. For research, such a tool could demonstrate the utility and impact of 
new vehicle technologies and guidance systems. 

1.1.1 Operational Applications 

In September of 2009 the Air Force’s new command, Air Force Global Strike 
Command, began operations to support the Global Strike mission. This command has 
taken over missions for the B-2 Spirit, B-52 Stratofortress and Minuteman III (LGM-30) 
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[1]. The mission of Global Strike is to be able to hit any desirable target in an 
advantageous time. The exact criterion is still flexible as it is more dictated by our vehicle 
limitations than flight planning limits. Though the Global Strike Command has specific 
units attached to its control, the concepts behind Global Strike are Air Force wide. The 
goal is to hit the enemy faster and harder from further away 

An easy to use mission planning tool would support this mission by allowing for 
new trajectories to be created rapidly. This would help the mission directors in planning 
attacks quickly and also permit them to make unlimited adjustments to meet demands. 
Being able to see numerous trajectory possibilities allows leaders to make well informed 
decisions quickly and match the appropriate weapon to the appropriate target with 
confidence. 

1.1.2 Air Force Research Lab’s Application 

The Air Force Research Laboratory (AFRL) is broken into groups called Focused 
Long-Term Challenges (FLTC). FLTC 4 is dedicated to Persistent and Responsive 
Precision Engagement, and within a subset of FLTC 4 lies challenge 4.2, Globally 
Deliver Full Spectrum of Kinetic Energy. The short-term goal for this challenge is to hit 
a target from anywhere in theatre subsonically. The mid-range goal is to hit a target 600 
nm in 20 min or 2000 nm in less than 2 hours at speeds greater than or equal to Mach 2. 

In the long-term, AFRL would like to be able to hit a target 9000 nm away in a time 
period of 2 hours or less. Eventually they would like to incorporate hypersonic speed and 
time periods of 1 hr or less to hit the desired target [8]. 
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As envisioned, a mission planning tool should help to facilitate the afore 
mentioned AFRL goals through its intended ability to set constraints and leave other 
variables free. If reaching a desired range is the goal other constraints can be allowed to 
change to see what is required to achieve that range. The same can be done with the time, 
speed and error. The speed at which trajectories can be calculated would allow for a 
greater depth of investigation in the same amount of time. 

1.2 Problem Statement 

The concept of a tool that can map trajectories seems simple enough, yet in reality 
the dynamics and constraints pose an intense problem to solve. This problem involves 
non-linear differential equations, two point boundary value problems, and path 
constraints. The top level problem description is to design an optimized trajectory to get a 
missile from a launch location to a prescribed target location. 

In order to accomplish the task of creating the optimized trajectory many subtasks 
need to be accomplished. Below are a few goals to solving this problem: 

1. Break the problem up into phases 

2. Define the reference frames and derive the equations of motion for each phase 

3. Connect the phases 

4. Set up constraints for the problem 

5. Develop a specific case to test and support the above goals. 

In order to solve the complicated problem of optimizing a trajectory, the optimal 
control problem needs to be broken up into 5 phases: 3 at launch (due to 3 stage vehicle), 
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a re-entry phase and a terminal phase. These phases are then optimized to minimize 
overall time in support of the Global Strike mission. 

After defining the phases, the reference frame for each phase has to be defined to 
allow for the phases to be tied together. The first section of phases (launch) is defined in 
an earth centered inertial reference frame and the other two phases are both in earth 
centered relative reference with the re-entry states defined in spherical coordinates, and 
the terminal phase in x, y, and z coordinates. 

The equations of motion for each section are based on the reference frame that the 
section is in. The first section has thrust, drag and gravity and has thrust control. The 
second section equations are the re-entry equations that include drag and lift and has bank 
angle as a control. The final section is assumed to be ballistic and has no control and is 
driven by gravity and the initial conditions for that phase. 

After all the phases are defined by their equations and constraints, the complete 
problem is given to the solver. The solver being used is General Pseudospectral 
Optimization control Software (GPOPS). This solver provides a solution for all of the 
above given information. This solver will output a trajectory, the control and Hamiltonian 
for the optimal control problem. 

1.3 Preview 

The rest of this thesis is broken down into four chapters; Chapter II will provide a 

literature review and background for the research. It will include a history of dynamic 

optimization and psuedospectral methods. It also includes information needed to design a 

mission planning tool such as: possible vehicles, atmospheric models and launch or target 
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sites. Finally, it includes a brief description of previous tools used to accomplish this goal 
and research that has been done in this area. Chapter III is the methodology and it covers 
the problem statement, assumptions crucial to the mission, how the problem is broken up 
and the method in which it is to be solved. Also included in this chapter is the break 
down as to how the mission planning tool will work when all the pieces are integrated 
correctly. Chapter IV is the results and analysis of the research. It includes the final 
trajectory and graphics which illustrate how different inputs impacted the trajectory. 
Finally in Chapter V contains the conclusions and suggestions for future research. 
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II. Literature Review 


2.1 Dynamic Optimization 

The concepts of dynamic optimization date back to the 17 th century when Johann 
Bernoulli posed a problem to his fellow scholars, such as Leibniz, THopital, Newton and 
Bernoulli’s brother, Daniel. This problem was called the brachistochrone problem which 
is Greek for ‘shortest time’. This problem was very true to its name, the goal was to find 
the minimum time for a bead to move down a string. From the idea of the 
brachistochrone problem a new field of mathematics was bom, optimization theory. In 
1733 Euler wrote Elmenta Calculi Variatonum, translated meaning Calculus of 
Variations and thus the founding mathematics behind optimization were established [12]. 

In the 1950’s Bellman made a contribution to Dynamic programming in updating 
the work of Hamilton and Jacobi thus fonnulating the Hamiliton-Jacobi-Bellman 
equation. In 1962 Pontryagin’s max(min) principle revolutionalized optimization theory 
in providing a method to detennine optimal control in a constrained problem. It is the 
underlying principle behind a bang-bang optimal control solution. 

Equation 1 shows the generic problem statement for a current day, optimal control 
problem, breaking it up into a cost equation, dynamics, boundary conditions and 
constraints. 
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rtf 

J = (p{x f ,t f )+\ L(x,u,t)dt 

Jt 0 

subject to the system dynamics: 
x ~-i(x,u,t) 

the boundary conditions: (1) 

t 0 and x 0 
i/s(x f ,tf) = 0 

and constraints (equality and inequality): 

C(x f ,t f ,x 0 ,t 0 ) < 0 

where J is the cost function, L is the Lagrange, and <p is the Mayer, x is the state vector, 
x is the dynamic vector (derivative of the states), t is the time, to and xo are the initial 
conditions, ip contains the function for the final conditions and C is the constraints. 

Analytical solutions to optimal control problems exist for simple cases. The 
majority of real-world problems would be too complex to solve using analytical methods 
and a numerical approach is needed. There are a numerous numerical methods in use 
today. 

Numerical methods can be either indirect or direct. Indirect methods use the first 
order optimality condition and solve the Hamiltonian for a boundary value problem. To 
solve problems indirectly, Calculus of Variations or the Method of Lagrange Multipliers 
are used. Two Point Boundary Value Problems can be solved with multiple standard 
methods. Small radii of convergence is a method that integrates states forward and the 
co-states backwards to solve for a solution. Another method, Steepest Decent, looks at 
the gradients of the dynamics and performance indexes to solve the differential equations 
[25]. 
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Direct methods are based on parameterizing functions with approximations and 
then transcribing them to a nonlinear program. These methods are the basis of 
pseudospectral methods. 

2.2 Pseudospectral Methods 

Pseudospectral methods are a subset of numerical methods. In order to explain 
Pseudospectral solution methods one must first understand nonlinear programming 
(NLP). NLP is the mathematical process of solving a system of equations incorporating 
both equalities and inequalities with a cost or objective function that is to be minimized 
or maximized. The nonlinear part is in the constraint equations or cost function. In 
equation fonn an NLP looks similar to the generic case below: 

Min f(z) (2) 

s.t. g(z)<0 
h(z)= 0 
zGZ 

where g is inequality constraint, h is the equality constraint, f is the objective/cost 
function, and z is a vector of parameters that exist in the domain of Z to include a set of 
solutions which satisfy all constraints. 

There are three different NLP solvers used regularly with MATLAB®: SNOPT, 

NPSOL and IPOPT. SNOPT “employs a sparse SQP algorithm with limited-memory 

quasi-Newton approximation to the Hessian or Lagrangian” whereas NPSOL employs a 

dense SQP algorithm. [2] Both SNOPT and NPSOL are effective on large scale 
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problems, but they have computationally expensive functions and gradients. IPOPT is 
also a large scale solver, but uses an interior point method. All of the NLP solvers use a 
nodal method to solve the problem. The problem in Equation 2 is broken up into phases 
and then discretized to a finite number of points (nodes) assigned to each phase and then 
the problem is solved at those points. An example of a multiple-phase problem is shown 
in Figure 1. 



I- 1-1 

Phase I Phase 2 Phase 3 Phase 4 


► 

lime 


Figure 1: Linkage for Multiple-Phase Optimal Control Problem [34] 


Next, a collocation scheme is used to capture the dynamics at each of the nodes. 
The approximations at each point come from 1 of 3 collocation schemes, Lobatto, Radau, 
or Gauss. A Lobatto scheme enforces conditions on the endpoint of the area of interest. 
An example of this method would be the trapezoid rule or Flermite-Simpson method. A 
Radau scheme enforces conditions only at one end of the area of interest. An example of 
this would be the Runge-Kutta method. Finally a Gauss scheme has conditions enforced 
strictly on the interior of the interval. An example of this method is the midpoint rule. 


9 






When a collocation scheme is combined with an NLP solver, a tool can be developed to 
solve nonlinear optimization problems. [33] 

There are multiple software packages that are used to solve optimal control 
problems. Built into MATLAB®, is a toolbox, ‘optimtooT which will solve both 
constrained and unconstrained problems as well as controlled and uncontrolled problems. 
However there are limitations to this software and it would be hard to set up a multiple 
phase problem. 

Software packages have implemented the pseudospectral method for solving 
optimal control problems. These are programs such as DIDO and GPOPS (General 
Pseudospectral OPtimal Control Software). Both of these will solve nonlinear optimal 
control problems with path constraints. GPOPS is a Gaussian Pseudospectral Method 
(GPM) program written in MATLAB® and solves complicated multiple phase optimal 
control problems. It’s structured in multiple levels to allow for each phase to have user 
defined specifications. GPOPS resembles other software in its structure that were written 
in FORTRAN as opposed to MATLAB®. 

2.3 Environmental Models 

In setting up any problem, including any of the methods explained in the previous 
section, the conditions need to be defined in which the problem is solved. Thus the 
environmental models used in this method need to be explored. The atmospheric and 
earth models are described below. 
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2.3.1 Atmospheric Models 

A basic atmospheric model would be a simple exponential model. This would 
have density and pressure as an exponential function of altitude. The most commonly 
used atmospheric model however is the standard atmospheric model. 

The standard atmosphere model is based off the hydrostatic equation. Assuming 
that the velocity is zero (no winds) a pressure, temperature and density model are 
developed as a function of altitude. The Mean Sea Level (MSL) values are from data 
collections: 

T 0 = 288.16 K = 518.69 °R 
P 0 = 101.325 kPa = 2116.2 lb/ft 2 
p 0 = 1.225 kg/m 3 — .002377 slug/ft 3 

The temperature model can be seen in Figure 2. It is approximated into areas of 
constant temperature (isothermal layers) and areas of varying temperature (gradient 
layers) based off of scientific data. 
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120 



Temperature (K) 


Figure 2: Temperature Profile for Standard Atmosphere [6] 

The equations for pressure and density are based off of the layer they fall in since 
they are reliant on the temperature. The equations for isothermal layers are equal because 
temperature is constant. Equation 3 shows the isothermal layers equation for density and 
pressure varying over altitude for a constant temperature. 

= L = e -( g / R T>- h i) (3) 

Pi Pi 

where h is the current altitude and hi is the reference altitude for the isothermal layer, p is 
the density at the current altitude and p i is the reference density for the isothermal layer, 
P is the pressure at the current altitude and P i is the reference pressure for the isothermal 
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layer, R is the universal gas constant, g is the gravity acceleration at the current altitude 
and T is the current temperature for the layer. 

For the gradient layers however temperature, pressure and density all have to be 
defined for the given altitude. Equations 4, 5, and 6 show the variance over altitude for 
the gradient layers. 


T^ + ThCh-hi) 

P _ /T\ _ ( g/ T h R) 

T 1 ~ vv 


_p_ 

pi 



t( g /T h R) + 1 ] 


( 4 ) 

( 5 ) 

( 6 ) 


with all the terms defined as in Equation 3 but T is the temperature at the current altitude 
and Ti is the reference temperature for the gradient layer, Th is the temperature gradient. 

For convenience the standard atmosphere model has been calculated and placed 
into tables. These tables have temperature, density, pressure and other properties verses 
altitude. The common used standard was published in 1976 which is contained in a table 
in the Appendix. Originally published in 1958 by the U.S. Committee on Extension to 
the Standard Atmosphere and updated in 1962, 1966 and finally 1976. Figure 3 shows 
temperature, density and pressure verses altitude plots. 
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The 1976 U.S. Standard Atmosphere 



temperature (K) 


density (kg/m 3 ) 


pressure (kPa) 


Figure 3: 1976 Standard Atmosphere Profiles [36] 


2.3.2 Gravity and Earth Model 

The two main models for the earth are a spherical one or an oblate spheroid one. 
Each of these models has a gravity model associated with it. The simplest model is that 
the earth is spherical. The spherical model means that the earth has a constant radius and 
has no flattening making the surface unifonn. This means that gravity can be measured as 
a function of radius from the center of the earth. 


E = - (7> 

8 r 2 

where p is the gravitational constant of earth and r is the radius from the center of the 
earth to the vehicle. 

This model assumes that the mass of the vehicle is negligible in comparison to the 
mass of the earth and gives a uniform density across a set altitude regardless of placement 
on the earth. 
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The other model having an elliptical surface was most recently defined by WGS- 
84 (World Geodetic System - 1984) [26], which is the currently valid scheme used by the 
Global Positioning System (GPS). This model has an oblate earth meaning that the poles 
are flattened. It has a major and minor axis for the earth which is placed into the model. 
MATLAB® currently has a built in function to use this model. The gravitation model 
associated with WGS 84 is the 1996 Earth Gravitational Model (EGM96). This model is 
designed as a group effort by the National Imagery and Mapping Agency, NASA 
Goddard Space Flight Center and the Ohio State University. Below is an image showing 
the geoid heights across the surface of the earth. 




EGM96 geoid heights in m. 


Figure 4: EGM96 Gravity Model [26] 


2.4 Highlighted Vehicles 

Along with multiple environmental models there are also multiple vehicles that 

could be used in this problem. These vehicles are explored below. 
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Figure 5: Minotaur III Launch [27] 


The Minotaur is an US rocket transformed from the LGM-118 Peacekeeper 
Missile. Its mission is to carry heavy payloads over a long range by a suborbital launch. 


Table 1: Minotaur III Data [28] 



Stage 1 

Stage 2 

Stage 3 


SR118 

SR119 

SR120 

Bum Time (sec) 

83 

54 

62 

Thrust (kN) 

1607 

1365 

329 

Mass Initial (kg) 

53070 

27800 

8200 

Mass Empty (kg) 

4086 

2900 

1400 

Diameter (m) 

2.36 

2.35 

2.35 

Length (m) 

7.72 

5.4 

2.3 

Shroud (stage 3) - 338 kg 

Payload - 1248 kg 

Nose Radius - 1 m 
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Figure 6: Space Shuttle Launch [35] 

NASA’s Space Shuttle went operational in the 1980’s and is now schedule to be 
retired in 2010. This vehicle’s missions include reaching an orbital trajectory and a re¬ 
entry process on its return, usually carrying people and supplies to and from the 
International Space Station. 

Though the Space Shuttle’s trajectory is not of great interest to the Air Force 
currently, the amount of data on its trajectories and heating profile is vast. This will allow 
for validations of the tool to be completed by making comparisons to the Shuttle model 
that will be discussed in Section 2.6 
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Table 2: Space Shuttle Data [35] 



Boosters 

External Tank 

Orbiter 


2 Boosters 

3 SSMEs 

2 OME 

Thrust 

2,800,000 lbf 
each 

1,225,704 lbf 
total sea level 

53.4 kN (12,0001bf) 
vaccum 

Specific Impulse 

269 s 

455 s 

316s 

Bum Time 

124 s 

480 s 

1250s 

Fuel 

Solid 

LOX/LH2 

MMH/N 2 O 4 

Height 

184 ft 

Diameter 

28.5 ft 

Mass 

4,470,000 lb 



Figure 7: Early artist’s rendition of X-37 [38] 

Designed by Boeing to test future launch technology for orbit and re-entry, the X- 
37 is a reusable robotic spacecraft originally designed to launch from the Shuttle cargo 
bay was redesigned for a Delta IV or comparable rocket. 


Table 3: X-37 Data [38] 


X-37 

Length 

27 ft 6 in 

Wingspan 

15 ft 

Weight (Loaded) 

12,000 lb 

Power 

Rocketdyne (Pratt & Whitney) 
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Figure 8: Simulated in-flight View of X-33 [37] 

Designed by Lockheed Martin as a tool to test capabilities required of a single 

stage to orbit reusable launch vehicle, the X-33 project was cancelled in 2001 but 

Lockheed Martin continues research on it. 


Table 4: X-33 Data [37] 


X-33 

Height 

20 m 

Mass 

285,000 lb 

Engine 

2 J-2S Linear Aerospikes 

Thrust 

410,000 lbf 

Fuel 

LOX/LH 2 
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Figure 9: Minuteman Launch [29] 

Currently the main missile for Global Strike Command, the Minuteman carries a 
nuclear warhead from a silo, using 3 stages to get to the target. 


Table 5: Minuteman Data [29] 



Stage 1 

Stage 2 

Stage 3 

Propulsion 

TU-122 

SR-19-AJ-1 

SR73-AJ/TC-1 

Weight 

78,0001b 

Length 

59 ft 9.5 in 

Diameter 

5 ft 6 in 

Range 

8100 miles 

Ceiling 

700 miles 

Speed 

15000 mph 
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2.5 Launch Sites 


Different vehicles are launched from different areas, thus it is important to 
examine diverse launch sites. It is important to have multiple launch sites in order to 
reach targets efficiently. Table 6 lists the most likely launch sites used by the United 
States. Figure 10 shows a map of all of the potential launch sites. The launch data came 
from previous cases ran by older tools. 



Figure 10: Map of Potential Launch Sites 


Table 6: Launch Site Data 


Launch 

Latitude 

Longitude 

Azimuth 

VAFB 

34.58 

-120.63 

170 

270 

KSC 

28.39 

-80.6 

35 

110 

WFF 

37.84 

284.51 

90 

160 

KLC 

57.44 

-152.34 

— 

— 

Equator 

0.00 

0.00 

170 

270 
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Vandenberg Air Force Base - VAFB 

VAFB is a military installation responsible for satellite launches for both military 
and commercial operations. It is also a test site for ICBMs (Intercontinental Ballistic 
Missiles) to include the Minuteman. 

Kennedy Space Center - KSC 

KSC is NASA’s main launch facility best known for housing and launching the 
Space Shuttle. Connected to KSC is Cape Canaveral Air Force Station. This has been and 
currently is the launch site for the Delta, Atlas and Titan rockets. In the future it will also 
be the launch site of SpaceX’s Falcon 9. 

Wallops Flight Facility - WFF 

WFF is located on the eastern shore of Virginia and is controlled by NASA 
Goddard as primary a research launch site, launching sounding rockets and small 
suborbital and orbital rockets. 

Kodiak Launch Complex - KLC 

KLC is located on Kodiak Island in Alaska and is controlled by the Alaska 
Aerospace Development Corporations. It is a commercial launch facility for suborbital 
and orbital launch vehicles. 

Equatorial Launch 

An equatorial launch site is ideal for putting satellites or orbiters into an 
equatorial orbit due to the fact that the vehicle can be directly launched into that orbit as 
opposed to having to perform a maneuver once in space and takes advantage of earth’s 
rotation minimizing fuel requirements. 
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2.6 Previous Models 


There are multiple methods that have been used in the past as well as current 


methods that generate trajectories for various vehicles. Some create optimized 


trajectories, others show feasible solutions. 


2.6.1 Shuttle Model 


The Space Shuttle model described below is applicable to this research because a lot 
of previous research for lifting bodies comes from Space Shuttle data. In order to 
compare the data the model must be understood. 

The trajectory is calculated from the drag acceleration envelope seen in Figure 11, 
which included limits such as: surface temperature, dynamic pressure, load, and glide 
limitations. [16] 


Along with implementing constraints in this model the trajectory is broken up into 


phases of flight. The three phases are: Temperature control, Equilibrium glide and 


Constant drag. 


DRAG 

ACCELERATION 

FT/SEC2 



RELATIVE SPEED THOUSANDS OF FPS 

Figure 11: Shuttle Trajectory Envelopes [16] 


23 



2.6.2 X-33 andX-3 7 Model 


Five main methods for trajectory calculations were tested and compared [15]. 
Their algorithms include: Baseline guidance [14], Linear quadratic method [11], 
Numerical predictor-corrector method [39], Drag-energy 3 dimensional method [30] and 
Quasi-equilibrium glide method [34]. The best being Quasi-equilibrium glide which led 
to Evolved Acceleration Guidance Logic for Entry (EAGLE) is shown in Figure 12. 



Figure 12: EAGLE Ground Track [33] 


Figure 13 shows how some of the previous models can use waypoint guidance to 
calculate the trajectory. The problem is broken into sections and guided from one section 
to the next solved not as a continuous problem but as a segmented optimization. 



Figure 13: Waypoint Guidance [30] 


24 





2.6.3 Currently Used Methods 

There are currently two methods used by the government to optimize trajectories, 
POST (Program to Optimize Simulated Trajectories) and OTIS (Optimal Trajectories by 
Implicit Simulation). POST was written by Lockheed Martin Astronautics in a joint 
effort with NASA Langley Research Center. POST is a two point boundary value solver 
that uses a direct shooting method to solve the equations of motion for the trajectory. This 
means that the solver takes a set of initial conditions and integrates forward in time to see 
if the final conditions are met. If they aren’t, then the initial conditions are changed and 
the integration is done again. This method is very computationally expensive thus taking 
a long time to find an optimized trajectory. POST is also not easily capable of 
incorporating waypoints or no-fly zones into trajectories. 

OTIS was written in the late 1980’s by Boeing and NASA Glenn Research 
Center. It has two methods for solving trajectories. It can do a direct shooting method 
similar to that of POST where it integrates forward in time from initial conditions. It also 
has a collocation method, where it fits a polynomial to the problem, allowing the 
trajectory to be solved at nodes along the trajectory as opposed to just at the endpoints. 

2.7GPOPS 

The tool used for this research is GPOPS (General Pseudospectral OPtimal 

control Software). This tool uses a nonlinear problem solver (SNOPT or IPOPT) to solve 

a set of differential equations. It has a collocation method allowing the problem to be 

solved at a user specified number of nodes. Though not specifically designed to solve 

trajectories, it is designed to solve multiple phase optimal control problems. The tool has 
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the user define a cost function, the dynamics or differential equations, connections 
between phases, and sets up limits and guesses for each of the states and controls. [4] 

Each GPOPS case is setup as a group of scripts. The four foundational scripts are: 
a main script; a DAE script (dynamics); a connect script; and a cost script. A brief 
description of each script is given below since these are key elements to the mission 
planning tool. 

Main Script 

The main script holds the framework for the case. It is setup into sections based 
off the number of phases a problem contains. A phase is defined by a distinction in the 
states, controls, equations of motion or cost functions. This structure allows the user to 
implement different phases to be used in one optimal control problem. Each phase has a 
select number of nodes assigned to it. Nodes are the number of points that will be passed 
into the NLP solver. This number is crucial because it must be large enough to catch the 
general solution but not too large to cause numerical divergence or have unneeded 
calculations. 

Each phase has a definition section in which minimums and maximums are 
defined for time, states, and controls. Both time and control allow for a minimum and 
maximum value for both the initial and final value. The states had minimum and 
maximum inputs for the beginning, along the section, and the end. A path constraint will 
be defined separately in the DAE script (described next), but the minimum and maximum 
values are set in the main script. 
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In order to solve the solution to the dynamics that will be defined in DAE, an 
educated guess is needed to input as a rough solution. This guess includes time, states and 
controls. It consists of as many points as the user has. The guesses need to be within the 
bounds previously set and their validity highly impacts the run time of the code to find a 
solution. Between each phase a connection matrix is defined which is passed into the 
Connection script and connects the states and controls of one phase to the next phase. The 
equations for this connection are in the Connection script. The Main script just sets the 
minimum and maximum values for the equations. 

DAE (Differential algebraic equations) Script 

The DAE script consists of the dynamics for the optimal control problem. Each 
phase can have its own dynamics or all of the phases could be under one set of 
differential equations. Each state requires a first-order differential equation to be provided 
to the NLP solver. These equations of motion can be very complicated or very simple and 
GPOPS can handle both. Along with the dynamics the DAE script includes any path 
constraints. There isn’t a limit to the number of constraints. These constraints can be 
written in terms of both states and controls. 

Connect Script 

The Connect script contains the final states from the phase before and the initial 
states for the phase after. From the before and after states, equations are written to 
transfonn one state to another or to allow for gaps in the stages between phases. The 
easiest example of this is in a multiple stage launch problem. At the end of each phase a 
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portion of the rocket is dropped and thus the mass will have a difference between the end 
of stage 1 or phase 1 and stage 2 or phase 2. 

Cost Script 

The Cost script creates the cost function for the optimal control problem. It 
includes 2 terms: the Lagrange and Mayer functions as defined in Equation 1. There can 
be a separate cost function for each phase, or a sum of weighted cost functions for each 
phase to create a total cost function. 

2.8 Previous AFIT Research 

Major Timothy Jorris completed a dissertation entitled “Common Aero vehicle 
Autonomous Re-entry Trajectory Optimization Satisfying Waypoint and No-Fly Zone 
Constraints” [22], His research included trajectory generation, waypoint satisfaction, 
threat or no-flyzone avoidance and used a GPM. William Karasz continued Major Jorris’ 
research with “Optimal Re-entry Trajectory Terminal State Due to Variation in Waypoint 
Locations”. [24] This research will be an expansion upon their work. Using Karasz’s 
non-dimensional equations of motion for the re-entry section allows all of his research on 
waypoints to be directly applied to that section of the mission planning tool. Their work 
on the application of GPOPS is also directly applied in this research. 

2.9 Summary 

Since the 17 th century optimal control problems have been posed in all shapes and 
sizes and this field has been shaped by famous mathematicians over the years. The 
advancements in this field have greatly been formed by the technological advancements 
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of the day. Starting with problems that can only be solved by hand to advanced nonlinear 
programming solvers, computerized tools are driving the current field of optimal control. 
The application of optimal control theory to trajectories seems natural. Optimizing 
trajectories ties with optimizing missile systems and a more optimized military. In 
civilian operations, optimized equates to less money. That’s why from the background 
research provided in the previous chapter, the mission planning tool objectives which will 
be described in the next chapter is a forward progression in optimal control. 


29 



III. Methodology 


3.1 Generic Problem Statement 

The overall objective of this research is to set up the foundations for a tool to 
compute the optimal trajectory. This foundation will calculate an optimized trajectory 
with settable parameters to allow for manipulation of the trajectory constraints. The 
trajectory will be broken up into phases with different equations of motion to get from a 
launch site to a target site. The goal is to obtain a rapid method for generating optimized 
trajectories. This will also show how the Gaussian psuedospectral solver, GPOPS, can be 
used on real life problems. 

This solver is set up in a modular setting in order to allow for multiple vehicles, 
launch and target sites, earth and atmospheric models and different constraints to be 
placed. The trajectory created will be optimized for minimum flight time. 

3.2 Mission Assumptions 

The assumptions for this mission are broken up into three categories: overall 
assumptions which apply to each phase for any configuration; phase assumptions that 
vary based on the phase the mission is in; and environmental assumptions. These are 
assumptions that are specific to an atmospheric or earth model. Listed below are the 
specific assumptions for each category. 

3.2.1 Overall Assumptions 

1. Vehicle is a point mass 

2. All mass separations are instantaneous 
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3. Lift and Drag Coefficients are constant 

4. Launch and Target Site Specified 

3.2.2 Phase Assumptions 
Launch 

1. Thrust is Constant - No throttling 

2. Only Control - 3 angles for thrust vectoring 

Re-entry 

1. No Thrust 

2. Rank angle only control 

Tenninal 

1. No Thrust 

2. No Control 

Model Assumptions 

1. There are only two atmospheric possibilities, the simple exponential and standard 
atmospheric model 

2. There are only two possible earth models, spherical and WSG-84 

3. The gravity models go with the earth model assumed (spherical to inversely 
proportional to radius squared and WSG 84 to EGM96) 
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3.3 Phases 



Launch 

3 Stages 
Inertial XYZ 


Re-entry 

Waypoints 
Relative R theta phi 


Terminal 

Single Phase 
Relative XYZ 


Minimize Time 


Figure 14: Phase Description 

The problem is broken up into 5 phases, 3 launch phases (1 for each stage), a re-entry 
phase and a terminal phase as shown in Figure 14. An exoatmospheric phase was considered, but 
initially to keep the problem simple, it is assumed re-entry conditions are met at the end of the 
launch phase. Thus, there are 3 sections in this problem launch, re-entry and tenninal. Each 
section has a different set of dynamics and is in a different reference frame. The launch section 
could be set up as 1 or 3 phases but 3 phases allow for easier transition between mass separations 
at each phase. There is one unifying cost function for the entire problem which is across all the 
phases, minimizing time. 
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3.3.1 Reference Frame 

There are four reference frames used in this mission. 

1. e is the inertial planet fixed reference frame 

2. e x is the relative planet fixed reference frame 

3. e 2 is the relative vehicle pointing reference frame 

4. e„ is the relative velocity reference frame 



Figure 15: Reference Frame Graphic [17] 


Reference frame 1, stays fixed in one spot regardless of time, thus the earth rotates 
around this reference frame. Reference frame 2 incorporates the rotation of the earth. This frame 
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stays fixed to the surface of the earth and is used to show the results. Reference frame 3 and 4 are 
show the angles for the vehicle and velocity directions. 

The phases are not all set in the same reference frame. The initial 3 launch phases are set 
in an inertial earth centered reference frame. The re-entry phase is in a relative earth centered but 
is in spherical coordinates and the terminal phase is also in relative earth centered coordinates 
but Cartesian x, y, and z. Below are the equations to transfonn each reference frame into a 
relative earth fixed frame (e,). 


ei 


1 0 
0 1 
0 0 


0 

0 

oo e At. 


( 8 ) 


e 


i — 


cos(0) cos (c(>) —sin (0) 

sin(0) cos (0) cos (0) 

sin (cj>) 0 


—cos(0)sin (4>) 
— sin(0) sin (4>) 
cos (0) 


e 2 


( 9 ) 



cos(0) cos (4>) 

—sin (0) 

—cos(0)sin (4>) 


cos (y) 

sin (y) 

0 

§1 = 

sin(0) cos (c()) 

cos (0) 

— sin(0) sin (40 


—cos(iJj) sin (y) 

COS (4^) COS (y) 

-sin (ijj) 


sin (4>) 

0 

cos (4>) 


— sin(if) sin(Y) 

sin(cp) cos(y) 

cos (40 


( 10 ) 


where 0, 0, y, and 0 are angles defined in Figure 15. <u e is the rotational rate of the earth and At 
is the time elapsed since e t equals e. 
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3.3.2 Equations of Motion 

One of the required MATLAB® scripts for GPOPS is DAE. DAE includes the dynamics 
of the problem which in this case would be the equations of motion (EOM). Each section has a 
different set of equations because they are in different reference frames and have different forces 
acting on them. Below are the 3 sections EOMs. 


Launch EOMs 

f — v (11) 

D + T (12) 

v =-g 

m 

m = constant (13) 

The states are: radius, r which consists of x, y and z, velocity, v which consists of v x v y and v z 
and mass, m. D is the drag, g is gravity and T is thrust. 

After the launch is completed the body has a constant mass, meaning mass is no longer a 
state in the equations of motion. In order to make the scaling for the problem easier for the solver 
to handle the dynamics are non-dimensionalized by the radius of earth, gravity at the surface or a 
combination of the two. The re-entry equations of motion have spherical state variables and are 
listed below. 
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Re-entry EOMs [24] 


r ’p = R V p siny 


r ■ gsiny 

V p = —D p -— + r p to e cos4>(cos4>sinY — sincj)sinv|;cosY) 

r p 


R VpCosYCosi]; 

p r p cos4> 


4> r 


R V p cosYsiru|; 


L p r p 

Y p = ——cosa - 


gcosy V r p 

_ —H-cosy + 2(n) e (cos4)cosi|j) + -y-a) e 2 cosc|)(cos4)COSY + sinct)sini|;sinY) 

Vpin Tp Vp Tp Vp 


= 


L p sina V p 


rco„ 


Vpcosy r p V p 


cosYcos\];tanc|) + 2to e (sin\|/cosc()tanY — sine])) — —-sincficoscficosij/ 

Vp co sy 


(14) 

(15) 

(16) 

(17) 

(18) 
(19) 


The p subscript indicates that the value is non-dimensional. The states are non- 
dimensional radius, r p and non-dimensional relative velocity, R V P . u> e is also non- 
dimensionalized. The remaining states are angles; longitude 0, latitude 0, heading angle \|/ and 
flight path angle y. D p is non-dimensional drag, L p is non-dimensional lift, w e is the rotation of 
the earth and the control is bank angle, o. 

After the re-entry phase the missile takes a ballistic dive to the surface for detonation. 
Those equations of motion follow. 
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Terminal EOMs 


r = v 

(20) 

D 

(21) 

^m“ g 


The states are: radius, r is a vector in the Cartesian coordinate frame consisting of x, y and z, and 
velocity, v which consists of v x v y and v z . D is the drag, m is the mass and g is gravity. 

3.3.3 Constraints 

In addition to the equations of motion, path constraints can be defined for each phase that 
the solution must satisfy. Maximum heating and acceleration are crucial to a mission planning 
tool as structures are only capable of handling certain temperature and structural loads. Below 
are the equations for heating and acceleration. 

<ts = n 1 / 2 T 3/2 (22) 


where q s is the stagnation heating, T equals kinetic energy and 77 equals non-dimensionalize 
altitude and are defined below. 

t _1/V\ (23) 

2 \go r o/ 


where V is velocity, g 0 is the reference gravity and r 0 is the reference radius. 

= pSCo (24) 

^ 2mp 

where p is the density, S is the reference area, C D is the coefficient of drag, m is the mass and 
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(3 is the scale height. 

The acceleration equations also use both the kinetic energy and the scale height listed in 
Equations 23 and 24. The acceleration is important for both structural loading of the missile itself 
and sensitivity of the payload. 


adecei = 2pg 0 r 0 riT 



2 


(25) 


where /? is the scale height, go is the reference gravity, ro is the reference radius, T is the kinetic 
energy, 17 is the altitude, C D is the coefficient of drag and C L is the coefficient of lift. [17] 


3.4 Optimal Control Problem 

The total optimal control problem ties the cost function with the equations of motion to 
get a total problem statement. GPOPS takes the cost function, dynamics and constraints and 
solves the optimal control problem. The cost function for the problem of a prompt strike is to 
minimize time, thus the complete optimal control problem can be written as follows: 

Min J = tf (26) 

s. t. 

x = f(x(t),u(t),t) 

C(x(t 0 ),t 0 ,x(t f ),t f ) < 0 


where x represents the first order derivative of the state vector, x is the states, u is the control, t is 
time t 0 is the initial time and tf is the final time. C represents the path constraints. 
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3.5 Specific Case 

To test the capability of this tool, a simple specific case is selected. This hypothetical 
case is launched from Kennedy Space Center and the desired target is Timbuktu (3° Longitude, 
16.5° Latitude). A spherical earth model is used with an inverse square law gravity model and 
the atmosphere is a simple exponential model. This case is a fairly simplified case but it is a 
building block to show future capabilities. The specific inputs for each phase are discussed in 
Chapter IV. 
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3.6 Mission Planning Tool 


Choose Launch Site 


Choose Target Site 


Choose Vehicle 


Choose Atmospheric 
Model 


Choose Earth Model 



Choose Constraints 


Create New 
Constraint 














Figure 16: Mission Planning Tool Flowchart 

Above is the flowchart for the mission planning tool. It shows the options the user has in 
the models and allows for future options to easily be added in. This program consists of a 
frontend that has the model, launch, target and vehicle selection along with any constraint 
placements. The backend allows for the user to select the data output they want as figures or 
variables and then the core section is the GPOPS code that the frontend automatically fonnats. 


41 













Figure 17, Figure 18, and Figure 19 show examples of the Graphic User Interface 
(GUI). This GUI follows the algorithm shown in Figure 16. Figure 17 is the main screen 
which allows the user to select or create the launch site, target, vehicle and constraints. 



Figure 17: GUI for Inputs 

Figure 18 is one of the subset menus. It allows the user to create a new launch site and 
when run, sets up all of the launch inputs for the GPOPS. Figure 19 shows the backend of the 
GUI. This is where the user selects what outputs they would like to be displayed to the screen 
and a filename for the data file. 
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Figure 18: GUI for Launch Site Inputs 



Figure 19: GUI for Output Selection 
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3.7 Summary 

This chapter provided an outline for the methodology followed to solve the optimal 
control problem presented. Analyzing the dynamics to evaluate the cost function over the set of 
constraints will give an optimized solution. Tying all of the methods together to apply to the 
specific case will give the desired results. Finally, the format of a mission planning tool was 
discussed so the application of this specific case to this model can be seen. The next chapter will 
go into the results for the specific case described in Section 3.5. 
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IV. Analysis and Results 


4.1 Changes in Specific Case 

Originally this research was set up to have a unified cost function for connecting phases. 
This method ran into a lot of difficulties with finding a convergent solution. The procedure to run 
an optimal control problem in GPOPS with varying states between phases had yet to be 
examined (it is accomplished later in the results). This problem has both time and non- 
dimensionalize time and states in both spherical and Cartesian coordinates. There are also 
multiple reference frames to be converted between. There are numerous possible problems that 
could be incurred in connecting phases with different states, for example scaling between could 
be a big problem (which is why non-dimensionalize dynamic equations were used). 

The current method used with re-entry problems in GPOPS to get a converged solution is 
to vary state constraints to allow for a “close” solution to be found and then iterating closer to the 
actual constraint. This method is hard to implement on 5 phases, especially when the states vary. 
More about how this could be accomplished is discussed in Chapter V. 

Due to the complications in the unified method, a sectional method was implemented to 
get solutions. This sectional method breaks the problem into 3 separate problems: launch, re¬ 
entry and terminal. The re-entry section still has the same cost function and setup as previously 
mentioned. The launch section has a new cost function which minimizes the distance to target 
seen in Equation 27. The terminal phase is solved by stepping forward in time to the target 
solution from the final value of the re-entry phase. 
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2 2 

Min J (Of ^target) ”1" (4*f ^target) 


( 27 ) 


where Of is the longitude at the end of the launch phase and 0 targe t is the target’s longitude and 
c p as latitude. 

As can be seen in Figure 20 the 3 sections are similar to the phases with different 
problem statements. Each section is still connected to the section before and after it, this is just 
accomplished outside of the optimal control problem. 




Launch 

3 Stages 
Inertial XYZ 
Maximize Downrange 




Re-entry 

Waypoints 
Relative R theta phi 
Minimize Time 



Single Phase 
Relative XYZ 
Time Stepped ODE 


Figure 20: New Section Break Down 
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Each phase of the GPOPS code requires minimum and maximum values for the states 


and controls and a guess for them as well. The guess can be the end points or any points along 
the path. It could even be an entire trajectory. The tables to follow show those values for this 
problem. 


Table 7: Launch States (Section 1) 


State 

Min 

Max 

Guess 

X 

-12756 km 

12756 km 

916.3 km 

y 

-12756 km 

12756 km 

-5535.1 km 

z 

-12756 km 

12756 km 

3033.5 km 

Vx 

-10 km/s 

10 km/s 

0.4036 km/s 

Vy 

-10 km/s 

10 km/s 

0.0668 km/s 

V z 

-10 km/s 

10 km/s 

0 km/s 

111 

90589 kg 

1194 kg 

90589 kg 


The guess, min and max values were the same for all three stages of the launch. The 
guesses are the inertial, initial conditions, the x, y, and z are the position of the launch site and 
V x , V y and V z are the rotation of the earth seen in Table 7. 


Table 8: Re-entry States (Section 2) 


State 

Min 

Max 

Guess 

Guess 

r p 

1 

1.2 

1.0191 

1.02 

v P 

.01 

1 

0.6953 

.9 

0 

-1.3064 

.0524 

-1.3064 

0.0472 

0 

0.2880 

0.4853 

0.4853 

0.2906 

'P 

-6.2832 

6.2832 

-0.1304 

-0.4952 

T 

-1.5533 

1.5533 

0.3289 

0.0007 


The guess in Table 8 was generated from the launch and terminal conditions. The values 
for Re-entry are earth centered relative. The max and min were decided by previous experience 
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with GPOPS on what conditions allowed for converged solutions and stayed away from 
singularities in the equations of motion. 


Table 9: Terminal Initial Conditions (Section 3) 


State 

IC 

X 

6113.16km 

Y 

289.02 km 

Z 

1830.53 km 

v x 

.103 km/s 

Vy 

-.958 km/s 

V z 

-.492 km/s 


Since the tenninal phase is a forward propagation of the equations of motion, the only 
values needed to solve the problem are the initial conditions shown in Table 9. The values for the 
Terminal phase are also earth centered relative. 

After finding converged solutions for the different sections, the re-entry and tenninal 
section were later linked in one minimum time problem. Creating only 2 sections for this 
problem a launch section and a re-entry thru tennination section. These results would not have 
been possible without first finding a possible solution in the 3 section method. 
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4.2 Results of Simulation Scenarios 

The new method has been explained in the Chapter 4.1, this section will explore the 
results of that method. 

4.2.1 Launch Results 


Radius vs Time - Launch 



Figure 21: Radius vs. Time for Launch Section 


The majority of the energy expended in launch is used to counter the weight of the 
vehicle and gain altitude. Figure 21 shows the climb in the launch phase. It can be infered from 
Figure 22 and Figure 23 that the first two stages are primarly used to get the vehicle off the 
ground. The manuevering for the desired direction is not clear until the third stage. Figure 21 
also shows that most of the altitude is gained in the third stage which would also show that it has 
the most flexibility to be manuevered. 
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Launch - x 





Figure 22: Position for Launch Section 

Vx - Launch 



time(s) 


Figure 23: Velocity for Launch Section 
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The controls for the launch section are the x, y and z components of the thrust unit vector. 
Thus providing the thrust vectoring for each stage. Figure 24 shows that the majority of the 
thrust was in the x direction for the largest portion of the launch. This seems logical because 
most of the thrust is required to get the vehicle to altitude as compared to changing its direction. 
Figure 25 and Figure 26 show the y and z components of the thrust. The most drastic thrust 
vectoring changes occur in the 3 ld stage. 


80 


Control 1 (x component) vs. Time 
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Figure 24: Control 1 Launch 
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Control 2 (y component) vs. Time 



5_i_i_i_i_i_ 

80 90 100 110 120 130 140 


0.5 -1-1-1-1-1-r 

0 - 


130 140 150 160 170 180 190 200 

time(s) 


Stage 1 


Stage 2 


Stage 3 


Figure 25: Control 2 Launch 
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Control 3 (z component) vs. Time 
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Figure 26: Control 3 Launch 
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x -|g' 3 Hamiltonian vs. Time - Launch 



Figure 27: Hamiltonian for Launch Section 


Figure 27 shows the Hamiltonian for the launch section. The Hamiltonian for the launch 
is approximately zero for all three phases which means that the controls and states for this 
section are feasible. For a fixed final time the Hamiltonian should be zero. 
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Mass vs time - Launch 


x 10 4 



Figure 28: Mass Profile 


Figure 28 shows the mass profile for the launch section. The mass state is a set profile 
regardless of the control or initial and final conditions. 
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4.2.2 Re-entry Results with Bank Control 

With the launch completed the next section is re-entry. The complete re-entry trajectory 
is show in Figure 29. As can be seen, the re-entry section of the trajectory is the longest portion 
by both time and distance. Re-entry trajectories tend to be oscillatory, which can be shown in 
the Figure 30. The velocity also oscillates with the dives and climbs of the trajectory which can 
be seen in Figure 31. 

The states appear feasible and the solution is converged. When the control is examined 
later in this section along with the Hamiltonian the feasibility of re-entry trajectory is questioned. 



Figure 29: Globe Trajectory - Re-entry Section 
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Radius vs Time - Re -entry 



Figure 30: Radius for Re-entry Section 


Velocity - Re-entry 



Figure 31: Velocity for Re-entry Section 
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Theta and Phi vs time - Re-entry 




Figure 32: Heading Angles for Re-entry Section 
The heading angles shown in Figure 32 have a smooth profile. The velocity angles, 

shown in Figure 33, presents the velocity angles, showing that the velocity is constantly 

changing direction in order to meet the requirements for the trajectory. This is due to solving the 

trajectory at a set number of points and then interpolating between them. 


gamma and psi vs time - Re-entry 




Figure 33: Velocity Angles for Re-entry 
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Control Vs Time 



Figure 34: Control for Re-entry 

Figure 34 shows an oscillating control profile. Since the only control was bank angle in 
this case and there wasn’t a bank rate constraint, the bank angle was able to jump between 
positive and negative 180 degrees. Physically the vehicle cannot be at a bank angle of negative 
180 degrees one second and then positive 180 degrees the next. There has to be a reasonable 
transitional period. Figure 35 supports the conclusions drawn from the control. The Hamiltonian 
doesn’t appear stable thus showing that the solution converged upon isn’t feasible. This means 
that the control is solved at select nodes for given values but isn’t a continuous profile. 


58 





































10 


Hamiltonian - Re-entry 



Figure 35: Hamiltonian Re-entry 
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4.2.3 Re-entry Results with Bank Rate Control 

Since the Hamiltonian didn’t converge to -1 and the control, bank angle wasn’t physically 
feasible it seemed necessary to increase the order or derivative of the control. The results to 
follow have a control of bank rate instead of bank angle. Figure 36 shows the new trajectory on a 
globe. To follow is a comparison of the ha nk angle control re-entry phase data and a 245 node 
bank rate control data. 



Figure 36: Globe Result for Ra nk Rate Control 
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Radius vs Time - Re -entry 




Figure 37: Radius and Velocity Difference in Controls 


As can be seen in Figure 37, both the velocity and radius profiles for a bank rate control 
have milder oscillations than that with only bank angle as the control. A converged solution 
could not be found for the bank angle control at as high of a nodal count as the bank rate solution 
which may account for some of the differences but the more stable profile appearance is likely 
due to the added degree of control. 
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Bank Angle Vs Time 



Figure 38: Bank Angle Differences 

Figure 38 shows how different having a rate as a control is to having the angle itself as 
the control. The ha nk angle appears to be more sinusoidal and the bank rate is more bang bang 
control as opposed to a bang bang control being the bank angle. The profile for the ha nk rate 
control is physically feasible for the missile to fly, thus, making this total trajectory more 
realistic. 
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Figure 39: Difference in Hamiltonian 

Figure 39 shows a Hamiltonian that is a lot closer to -1 for the bank rate control than the 
bank angle control. A close look at the Hamiltonian for the bank angle rate is shown in Figure 


40. 
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Hamiltonian - Re-entry 



Figure 40: Hamiltonian for Bank Rate 

Control Vs Time 



Figure 41: Bank Rate vs. Time 
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The data compared in the comparison graphs for the ba nk rate control verses bank angle 
control are for 245 nodes. The number of nodes a problem uses drastically impacts the solution. 
A converged solution is identified by the same solution being presented for a various number of 
nodes. This leads me to believe that the Hamiltonian for the ha nk rate control could reach a 
more stable value as its convergence to -1 can be shown in Figure 46. Both Figure 42 and Figure 
43 show small difference in the different nodes the same is also true for Figure 44 the position. 
The bank angle has slightly greater differences because it is directly tied to the control bank rate. 
These results show a converging solution so one could say the re-entry phase is feasible. 


Radius vs Time - Re -entry 



Figure 42: Radius for Different Nodes 
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Velocity - Re-entry 



Figure 43: Velocity for Different Nodes 

Theta and Phi vs time - Re-entry 




Figure 44: Position for Different Nodes 
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Hamiltonian 


Bank Angle vs Time 



Figure 45: Bank Angle for Different Nodes 
Hamiltonian - Re-entry 



Figure 46: Hamiltonian for Different Nodes 
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4.2.4 Terminal Results 


The terminal results are purely created from a forward propagation of the initial 
conditions. Figure 47 shows a 3D representation of this section of the trajectory. 


3D plot for Terminal Section 



Figure 47: 3D position plot for Terminal Section 
The states are shown in Figure 48 and Figure 50 with Figure 49 showing the magnitude 

of the velocity. 
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Figure 48: Position for Terminal Section 
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Figure 49: Velocity Profile for Terminal Section 
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Figure 50: Velocity Breakdown for Terminal Section 
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4.2.5 Total Trajectory Results 



Figure 51: Globe Complete Trajectory 

Figure 51 shows the entire trajectory on earth from KSC to Timbuktu. Figure 52 and 
Figure 53 show the radius and velocity for the total trajectory showing smooth connection from 
one section to the next. There are no discontinuities in the plots. These plots use the bank rate 
control thus having a feasible solution. Yet it can be seen by the amount of maneuvering at the 
end of the trajectory that this solution is not optimal. 
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Radius vs Time for Total Trajectory 



time(s) 

Figure 52: Radius for Complete Trajectory 

Velocity vs Time for Total Trajectory 



Figure 53: Velocity for Complete Trajectory 
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4.2.6 Connected Results 


Due to the lack of optimality in the 3 separate solution method it seemed crucial to return 
to a connected approach. Instead of attacking the entire problem again, Only the re-entry phase 
and the terminal phase are connected. An initial converged solution took about 20hrs to achieve. 
The results create a clearly more optimal trajectory as can be seen in the figures to follow. Figure 
54 shows the connected radius and velocity for the re-entry and terminal phases. The radius and 
velocity both connect smoothly between phases and still trend towards similar patterns, but the 
total time is now around 1600 seconds for re-entry and terminal phase as opposed to over 2500 
seconds for just the re-entry. 


Radius and Velocity Vs Time 




Figure 54: Radius and Velocity for Connect Solution 
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Position Vs Time 




Figure 55: Longitude and Latitude for Connected Solution 


Velocity Angles Vs. Time 




Figure 56: Gamma and Psi for Connected Solution 
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Figure 57: Globe for Connected Solution 

Figure 55 shows a direct path in both longitude and latitude to the target site. This seems 
appropriate for a minimum time problem between two points to create a straight line; this can 
also be seen in the globe trajectory plot, Figure 57. Figure 56 shows a small variance in the 
velocity angles. 

Figure 58 shows the bank angle and ba nk rate for the re-entry portion of the connected 
solution. A saw tooth wave is the bank angle solution with a square wave as the control or ha nk 
angle. The other solutions in the re-entry phase appeared to be trending towards this solution but 
connecting the phases allowed for the solution to completely converge. Figure 59 shows the 
Hamiltonian is within .01 on either side of -2 as opposed to the pervious Hamiltonians which 
varied at least 10 times that. 
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Hamiltonian 


Bank Angle and Rate Vs Time for Re-entry 




Figure 58: Bank Angle and Rate for Re-entry Connected Solution 


Hamiltonian for Re-entry 



Figure 59: Hamiltonian for Connected Solution 
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4.3 Investigative Questions Answered 

Though the methodology did change in the process of this experiment, the original goals 
were still achieved. Below are a recap of the original goals and a brief description of how they 
were achieved. 

4.3.1 Break the problem into phases 

The problem was broken up into 3 main sections and then additional phases are added for 
stages in launch or waypoints in the re-entry section. Thus 5 total phases are assigned, for an 
efficient transition between different dynamics and states in the problem. 

4.3.2 Define the reference frames and derive the equations of motion for each phase 

The reference frames for each section are clearly defined. The process for converting 
from one reference frame to another is in place for any future cases. Each of the 3 sections has 
their own equation of motion defined for the reference frame the state will be in. 

4.3.3 Connect the phases 

The code to switch between reference frames has successfully been tested outside of the 
optimal control solver, showing that the method behind the connection equations is valid. The 
derivation for the connection equations can be found in Appendices A and B. 

4.3.4 Set up constraints for the problem 

All the constraints or fonnats for potential additional constraints have been set up for 
this problem. Constraints for the states and controls have been put into place by establishing the 
minimum and maximum values which are written in the Main script. The path constraints which 
would include acceleration and heating have not yet been introduced but the structure is there 
once minimums and maximums for the case are determined. 
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4.3.5 Develop specific case to test and support the above goals 

Section 4.2 provided a description of a trajectory for a specific case. All of the results are 
based upon this case. 

4.4 Summary 

Due to complication in the implementation of a single solution method, a sectional 
method was employed breaking the problem into 3 smaller problems: launch, re-entry, and 
tennination. The solutions for the section method are shown in Section 4.2. Due to infeasibility 
in the original setup, an additional degree of control was added to the problem, bank rate. These 
results then allowed for the re-entry and tenninal phases to be connected creating now only 2 
separate parts of the problem. Though not fully meeting the overall objective of designing and 
implementing a single solution method, all of the goals for setting up a mission planning tool 
were met, shown in Section 4.3. 
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V. Conclusions and Recommendations 


5.1 Conclusions of Research 

The results obtained demonstrate that to get a complete optimal trajectory for the enter 
mission as opposed to a piece-wise optimal solution for each section, the 3 separate trajectory 
solutions need to be tied together into a single cost and then solved by the GPOPS code. The 
current sing section trajectory solution can be used as a good guess to place into a complete 
GPOPS code at a later stage. The research show that for the re-entry section an addition state for 
a ba nk angle was needed for the control variable to keep a realistic bank angle. Overall, it was 
shown that GPOPS is capable of generating piece-wise missile trajectories and has potential for 
use as a mission planning tool. 

5.2 Significance of Research 

There is a substantial requirement for a new trajectory optimization tool in the Air Force. 
Computational abilities have drastically improved since 1988 when OTIS was contracted. The 
idea of rapid trajectory generation has changed over the years from days to minutes to the 
potential of seconds to create a trajectory. GPOPS has the potential to solve the complicated 
equations of motion faster than its predecessors thus making it ideal for the mission planning 
tool. This research is a great first step to a user friendly mission planning tool. It shows the 
mathematical capability of optimizing the stages of the problem in a rapid solution and setting up 
for future models to be tested. 
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5.3 Recommendations for Action 

The research provided in this paper leads to instant actions, or actions that can take place 
without much further research. The three main recommendations for action are: Integrate new 
GPOPS version 3.0; Compare case with POST or OTIS data; and run other cases of interest. 
GPOPS version 3.0 includes a new mesh grid algorithm which should no longer require the user 
to figure out the correct number of node points needed to solve the trajectory. This should allow 
for a more automatic process in solving trajectories. 

Since POST and OTIS are both currently being used today an immediate action to be 
perfonned would be to run similar scenarios on both GPOPS and POST or OTIS and compare 
the results and run time. GPOPS should be able to process an answer quicker but the differences 
between the solutions would need to be further investigated. 

This research was sponsored by AFRL with the specific case being a research case they 
have ran with their systems before, but simplified. Increasing the model to match AFRL’s cases 
and getting specific targets of significant to place in the model would make the trajectories easier 
to compare to historical data. In addition to actions that can currently be completed there are 
many additions to this research that would be worthy of future investigation. The next section 
covers those research endeavors. 

5.4 Recommendations for Future Research 

The main motivation behind this research is to create a mission planning tool so future 
research should be in support of that goal. There are many steps that are needed in order to get a 
user friendly mission planning tool that will rapidly generate trajectories. The future areas of 
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research for a mission planning tool can be broken up into five categories: Vehicle 
Specifications; Model; Path; Code and Far Future additions. 

Vehicle Specifications 

In regards to the aerodynamics information there is a lot of future work that can be done 
in this area. First, implementing an angle of attack dependency for the lift and drag coefficients. 
Also, not assuming that the mass is a point mass and thus has a changing center of mass and 
pressure with appropriated lift and drag and moment calculations. Since a lot of aerodynamic 
data is experimental and thus tabular, having a table lookup option for aerodynamics would be a 
nice addition to the vehicle models. 

In the propulsion section of the vehicle model, having an option for varying types of 
propulsion systems would be a needed expansion. Allowing for throttled thrust, variable mass 
flow rates or a specific thrust profile to be input would create more realistic and versatile models. 

The controls in the original set up are fairly simplistic having only angles or angular rates 
as the control. Creating options for more complicated controls and feedback loops into the 
dynamics that allow for reactions to onboard sensors or GPS data would create a more realistic 
solution. Also having a library of vehicles and data available for the tool would be very 
appealing to users to be able to select a vehicle as opposed to finding all the data for any new 
case one wanted to run. 

Model 

In addition to the atmospheric and gravity models, having wind models would increase 
the accuracy of the final trajectory. Though a lot of tools neglect wind to simplify their models, 
there are multiple standard wind models that could be apply to this tool to be selected by the 
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user. Also the number of models for both the atmosphere and gravity could be increase to offer 
the user of the tool more possibilities for research. One example would be to add a nonstandard 
atmospheric model. 

An interesting study that could possibly help make the overall tool run faster would be to 
analyze the speed of convergence for each additional level of fidelity in the model. Future 
research could determine whether or not running a case with simplifier models first to determine 
a guess solution would be faster than trying to get the complicated models to converge from a 
basic guess. 

Path 

Studying No-fly zones and Waypoints have been the primary research emphasis for many 
of my predecessors in this field. Applying what they have done for a mission planning tool and 
expanding on it would create highly developed path constraints. An ideal user function for the 
mission planning tool would be to have the user capable of selecting a treaty they would like to 
follow which would include a set of no-fly zones. 

Another path changing constraint that should be added to this tool is angle requirements 
for both launch and detonation. A target may need to be attacked from a specific direction in 
order to minimize civilian causalities, for example to avoid a hospital the missile made need to 
attack from a direction other than the optimized one. As could be seen with the launch site data, 
most launch sites have angle restrictions on the azimuth angle. Future research could be to 
detennine the impact changing these angles has on the entire trajectory. 
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Code 


The areas of coding that would benefit from future research would be to unify the 
sections created into a single case like originally planned. Preliminary steps would need to be 
taken to get a single code to work with different states. There may be scaling issues with the 
original set up of this problem. Future investigations may need to figure out the best method to 
create a single code with a unified cost function. 

A guess generator is also a great need for the mission planning tool. In order for the code 
to run the fastest, which speed is one of GPOPS greatest selling points for this tool, a better guess 
is needed than initial and final conditions. This guess generator could use an ordinary differential 
equation solver to gain a rough trajectory to be input into the guess for GPOPS. 

The current phase connections equations have to call multiple functions and sometimes 
calculate values that aren’t needed as a check for validity. A more efficient phase connection 
may also help to connect the entire code. Also increasing the control in the re-entry section to 
bank acceleration should provide smoother bank angle plots and may make the solution converge 
at a faster rate. 

Far Future 

Along with finding the optimal trajectory, it is important to have infonnation on how 
variations in the parameters will cause the mission to fail. The capability to do failure analysis 
would be a great future addition to the mission planning tool. Keeping the final conditions to 
within a set error tolerance minimum and maximum variations at waypoints or at launch could 
be calculated. Thus the tool would be creating an optimal trajectory with a success envelope 
around it. 
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This tool should be versatile enough to be adapted for many missions in the future and 
not just a missile trajectory. Orbital missions should not be too complicated to add on the 
mission. In addition to the transfonnation between reference frames that is already in place, 
orbital parameters would need to be added to the mix. 
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Appendix A: Rpg2Xyz 


The derivation below shows how position is converted into an earth fixed relative frame 


from a vehicle fixed frame. 
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Appendix B: Vgp2Vxyz 


The derivation below shows how velocity is converted into an earth fixed relative frame 


from the velocity frame. 
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Appendix C: 1976 Standard Atmosphere Table Example [33] 
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Appendix D: Code 


Example of Main Code 


o, _ 

o 

% Main File 
% Re-entry Section 

O, _ 

o 



% Danielle Yaple 
% 28 December 2009 


o, 

o 


function [Re entry] =Reentry Main() 

global CONSTANTS 

a = 100; 

nodenum = a; 

run MOD_basic 

run VEH minotaur 

run TAR_timbuktu 

run LAU_KSC 

CONSTANTS.m = 1194.8; %kg 
% rO = 6.5000e+003; 

% VOmax = 5.5; 

% thetaO = -1.3685; 

% phiO = 0.4950; 

% gammaO = -.6632; 

% psiO = -0.0279; 

rO = 6.500e+003; 

V0 = 5.5; 

thetaO = -1.306357339705933; 

phiO = 0.485335659325143; 

gammaO = 0.328910761317162; 

psiO = -0.130375293982734; 

sigma_guess = 0; 
gs = 9.81/1000; %km/s A 2 
rs = CONSTANTS.rs; 

CONSTANTS.wp earth = 7.2722E-5*sqrt(rs/gs); % rad/s 
rpO = rO/rs; 

VpO = VO/sqrt(rs*gs); 

thetaf = CONSTANTS.TARGET.theta; %rad 
phif = CONSTANTS.TARGET.phi; %rad 

%extremes 

tp_min = 500/sqrt(CONSTANTS.rs/gs); 
tp_max = 15000/sqrt(CONSTANTS.rs/gs); 
rpmin = 1; 
rpmax = 1.2; 
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Vpmin = 0.01; 

Vpmax = 1; 
thetamin = thetaO; 
thetamax = thetaf; 
phimin = phif; 
phimax = phiO; 
gammamin = -deg2rad(89); 
gammamax = deg2rad(89); 
psimin = -2*pi; 
psimax = 2*pi; 

sigmamin = deg2rad(-60); %rad 
sigmamax = deg2rad(60); %rad 


rp_guessl = 1.02; 
Vp_guessl = .9; 
gamma_guessl = 0; 
psi_guessl = deg2rad(0); 


%from ode 

rf = 6.387892378315309e+003; 
rfmax = rf; 

rpf = rf/CONSTANTS.rs; 
rpfmax = rfmax/CONSTANTS.rs; 
phif = 0.290637756866697; 

thetaf = 0.047243562179519; 

Vf = 1.082251805453242; 

Vfmax = Vf +1; 

Vpf = Vf/sqrt(CONSTANTS.rs*CONSTANTS.gs); 

Vpfmax = Vfmax/sqrt(CONSTANTS.rs*CONSTANTS.gs); 

% Vpfmin = Vfmin/sqrt(CONSTANTS.rs*CONSTANTS.gs) 
% Vpfmax = Vfmax/sqrt(CONSTANTS.rs*CONSTANTS.gs) 


gammaf 

= 7.051 

psif = 

-0.495 

iphase 

= 1; 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 

limits 

(iphase) 


.nodes = nodenum; 


.time.min = [0 
.time.max = [0 
.state.min(1, 

.state.max(1, 

.state.min(2, 

.state.max(2, 

.state.min(3, 

.state.max(3, 

.state.min(4, 

.state.max(4, 

.state.min(5, 

.state.max(5, 

.state.min(6, 

.state.max(6, 
.control.min(1, 
.control.max(1, 


tp_ 

tp_ 


maxj ; 

= [rpO rpmin rpf]; 

= [rpO rpmax rpf]; 

= [VpO Vpmin Vpf]; 

= [VpO Vpmax Vpf]; 

= [thetaO thetamin thetaf]; 

= [thetaO thetamax thetaf] ; 

= [phiO phimin phif]; 

= [phiO phimax phif]; 

= [gammamin gammamin gammaf] 
= [gammamax gammamax gammaf] 
= [psiO psimin psif]; 

= [psiO psimax psif]; 

= sigmamin; 

= sigmamax; 
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limits(iphase).parameter.min 
limits(iphase).parameter.max 
limits(iphase).path.min = []; 

limits(iphase).path.max = []; 

limits(iphase).event.min = []; 

limits(iphase).event.max = []; 


.mat' ], 'file' ) 

output.solution(iphase).time; 
output.solution(iphase).state; 

= output.solution(iphase).control; 

output.solution(iphase).parameter 


if exist([mfilename, 
load(mfilename) ; 
guess(iphase).time = 
guess(iphase).state = 
guess(iphase).control 
guess(iphase).parameter = 

else 

guess(iphase).time = [0; 

guess(iphase) .state (:,1) = 
guess(iphase) .state (:,2) = 
guess(iphase) .state (:,3) = 
guess(iphase) .state (:,4) = 
guess(iphase) .state (:,5) = 
guess(iphase).state(:,6) = 
guess(iphase) .control (:,1) 
guess(iphase).parameter = 

end 


tp_max]; 

[rpO;rpguessl]; 

[VpO;Vp_guessl]; 

[thetaO;thetaf]; 

[phiO; phif]; 

[gammaO; gamma guessl]; 

[psiO; psi_guessl]; 

= [sigma_guess; sigma_guess]; 


setup.name = 'Re-entry'; 

setup.funcs.cost = 'Reentry_Cost' ; 

setup.funcs.dae = 'Reentry_Dae' ; 

% setup.funcs.link = 'Reentry_Connect' ; 
setup.limits = limits; 
setup.guess = guess; 

% setup.linkages = linkages; 
setup.derivatives = 'automatic'; 
setup.direction = 'increasing'; 
setup.autoscale = 'on' ; 

output = gpops(setup); 
solution = output.solution; 
status = output.SNOPT^info; 

Re_entry.sol = solution; 
if status == 1 

save(mfilename, 'output' ) ; 
r = output.solution.state(:, 1)*rs ; 

V = output.solution.state(:,2)*sqrt(rs*gs) ; 

theta = output.solution.state(:, 3) ; 

phi = output.solution.state(:, 4) ; 

gamma = output.solution.state(:,5); 

psi = output.solution.state(:, 6) ; 

time = output.solution.time*sqrt(rs/gs); 

linestyles = ['m-']; 

plot_on_sphere(phi,theta,r,21,linestyles) ; 
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figure (22) ; 

subplot(7,1,1); plot(time,r-CONSTANTS.rs,linestyles) ; 
ylabel ({' altitude ','[ km] '}) ; axis tight; hold on; 

subplot(7,1,2); plot(time,V,linestyles); ylabel( { 'velocity' , ' [km/s] ' }) 
axis tight; hold on; 

subplot(7,1,3); plot(time,(theta).*180/pi,linestyles); 
ylabel ({' longitude ',' [deg] '}) ; axis tight; hold on; 

subplot(7,1,4); plot(time,(phi),*180/pi,linestyles); 
ylabel ({' latitude ',' [deg] '}) ; axis tight; hold on; 

subplot(7,1,5); plot(time,(gamma).*180/pi,linestyles); ylabel({ 'flight 
path-angle', '[deg] '}); axis tight; hold on; 

subplot(7,1,6); plot(time,(psi).*180/pi,linestyles); 
ylabel ({' heading ',' [deg] '}) ; axis tight; hold on; 

subplot(7,1,7); plot(theta.*180/pi,phi.*180/pi,linestyles); 
ylabel (' latitude [deg]'); xlabel (' longitude [deg]'); axis tight; hold on; 

1 = length(time); 
time(1); 
rf = r(1) ; 

Vf = V(1) ; 
thetaf = theta(1); 
phif = phi(1); 
gammaf = gamma(1); 
psif = psi(1); 

[x,y,z] = rpt2xyz(rf, phif, thetaf); 

[Vx,Vy,Vz] = Vgp2Vxyz(Vf,gammaf,psif,thetaf,phif); 


else 

fprintf( 'boo\n' ) 


end 
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Example of Connect Code 

function [connect Dconnect] = SimlConnect (sol); 

global CONSTANTS 

xf_left = sol.left.state; 
p left = sol.left.parameter; 
left_phase = sol.left.phase; 
xO_right = sol.right.state; 
p right = sol.right.parameter; 
right_phase = sol.right.phase; 

connect = xO_right-xf_left; 


% avoid calc of derivs in not necessary 
if nargout == 2 

Dconnect = [-eye(length(xf_left)), eye(length(xO_right))] 

end 
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Example DAE Code 


function XDOT = SimlDae(sol) 

% These are the equations of motion for re-entry based on Hicks 
% t = solode.time; 

rp = sol.state(:,1); % altitude ; rp = radius prime (prime to denote 
nondimensional) 

Vp = sol.state(:,2); % velocity-like ; Vp = velocity prime 

%theta = solode.states(:,3); % longitude 

phi = sol.state (:,4); % latitude 

gamma = sol.state(:, 5) ; % flight-path angle 

psi = sol.state(:,6); % heading 

% sigma = sol.state (:,7); 

% bnkrate = sol.control(:,1); % bank angle 
sigma = sol.control(:,1); 
ul = cos(sigma); 
u2 = sin(sigma); 

% constants ; the calling function should have a global variable named 
CONSTANTS with the below variables, 
global CONSTANTS 

rs = CONSTANTS.rs; % radius at the surface of the Earth [km] 

rhos = CONSTANTS.rhos; % density of atmosphere at surface of the Earth 

beta = CONSTANTS.beta; % scaling constant for atmosphere 

m = CONSTANTS.m; 

Cl = CONSTANTS.Cl; 

Cd = CONSTANTS.Cd; 

S = CONSTANTS.S; 

wp earth = CONSTANTS.wp earth; 

Lp = .5.*rhos.*exp(-beta.*rs.*(rp-1)).*C1.*S.*Vp. A 2.*rs./m; 

Dp = .5.*rhos.*exp(-beta.*rs.*(rp-1)).*Cd.*S.*Vp. A 2.*rs./m; 

drpdot = Vp.*sin(gamma); 
dVpdot = -sin(gamma)./rp. A 2 + 

rp.*wp earth. A 2.*cos(phi).*(cos(phi).*sin(gamma)- 
sin(phi).*sin(psi).*cos(gamma)) - Dp; 

dthetapdot = Vp./rp.*cos(gamma) .*cos (psi) ./cos(phi); 
dphipdot = Vp./rp.*cos(gamma).*sin(psi); 

dgammapdot = .5.*rhos.*exp(-beta.*rs.*(rp-1)).*C1.*S.*Vp.*rs.*ul./m - 
cos(gamma)./(rp. A 2.*Vp) + Vp.*cos(gamma)./rp + 

2.*wp earth.*cos(phi).*cos(psi) + 

rp./Vp.*wp earth. A 2.*cos(phi).*(cos(phi).*cos(gamma)+sin(phi),*sin(psi).*sin( 
gamma)); 

dpsipdot = Lp.*u2./Vp./cos(gamma) - Vp./rp.*cos(gamma).*cos(psi).*tan(phi) + 

2.*wp earth.*(sin(psi).*cos(phi).*tan(gamma)-sin(phi)) - 

rp.*wp earth. A 2./Vp.*sin(phi).*cos(phi),*cos(psi)./cos(gamma); 

XDOT = [drpdot dVpdot dthetapdot dphipdot dgammapdot dpsipdot]; 
end 
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Example of Cost Code 


function [Mayer,Lagrange]=Reentry_Cost(sol) 


tO = sol.initial.time; 
xO = sol.initial.state; 
tf = sol.terminal.time; 
xf = sol.terminal.state; 
t = sol.time; 
x = sol.state; 
u = sol.control; 

%Method 1 Minimize Time 

Mayer = tf; 

Lagrange = zeros(size(t)) ; 
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